# Probit Regressions for Response Variables using Pooled Cross-sectional Dataset

# LPG Connection (Binary)
model.dd.vfe.uselpg.scst_ols = fixest::feols(uselpg ~ scst*wave + age + gender +
                                     religion + education + hhsize +
                                     lnexp + bplaay | factor(village), 
                                   data = pool, weights = pool$weights)

model.dd.vfe.uselpg.scst_glm = glm(uselpg ~ scst*wave + age + gender +
                                           religion + education + hhsize +
                                           lnexp + bplaay + factor(village), 
                                         data = pool, weights = pool$weights,
                                         family = quasibinomial(link = "probit"))

coefs.dd.vfe.uselpg.scst_glm <- coeftest(model.dd.vfe.uselpg.scst_glm,
                                         function(x) sandwich::vcovHC(x))

# Grid Connection (Binary)
model.dd.vfe.usegrid.scst_ols= fixest::feols(usegrid ~ scst*wave + age + gender +
                                  religion + education + hhsize +
                                  lnexp + bplaay | factor(village), 
                                data = pool, weights = pool$weights)

model.dd.vfe.usegrid.scst_glm = glm(usegrid ~ scst*wave + age + gender +
                                            religion + education + hhsize +
                                            lnexp + bplaay + factor(village), 
                                          data = pool, weights = pool$weights,
                                          family = quasibinomial(link = "probit"))

coefs.dd.vfe.usegrid.scst_glm <- coeftest(model.dd.vfe.usegrid.scst_glm,
                                          function(x) sandwich::vcovHC(x))


models.vfe.dd_access_scst_ols <- list(model.dd.vfe.uselpg.scst_ols,
                                      model.dd.vfe.usegrid.scst_ols)

models.vfe.dd_access_scst_probit <- list(model.dd.vfe.uselpg.scst_glm,
                                          model.dd.vfe.usegrid.scst_glm)

coefs.vfe.dd_access_scst_probit <- list(coefs.dd.vfe.uselpg.scst_glm[,2],
                                        coefs.dd.vfe.usegrid.scst_glm[,2])

fixest::esttex(models.vfe.dd_access_scst_ols, 
       se = "cluster", 
       digits = 3,
       dict = c(scst = "SC/ST", wave = "2018", uselpg = "LPG Connection",
                usegrid = "Grid Connection", convenience = "LPG Home Delivery",
                `factor(village)` = "Village", hh = "Household"),
       keep = c("SC/ST", "2018"),
       fitstat = c("r2"),
       fixef_sizes = T,
       float = F,
       replace = T,
       file = c("./Manuscript/Tables/dd_access_scst_robustness_ols.tex"))

omit.stats <- c("adj.rsq","ser","f")

#Generate latex table
stargazer(se = coefs.vfe.dd_access_scst_probit,
          header = F,
          models.vfe.dd_access_scst_probit,
          keep = c("scst","wave","scst:wave"),
          covariate.labels = c("SC/ST", "2018", "SC/ST * 2018"),
          no.space = T,
          suppress.errors = F,
          float = F,
          omit.stat = omit.stats,
          dep.var.labels = c("LPG Connection",
                             "LPG Home Delivery",
                             "Grid Connection"),
          add.lines = list(c("Village FE", c("Yes","Yes","Yes"))),
          notes = c("1. *p<0.1; **p<0.05; ***p<0.01"),
          notes.append = FALSE,
          notes.align = "l",
          out = c("./Manuscript/Tables/dd_access_scst_robustness_probit.tex"))

margins_uselpg <- margins(model.dd.vfe.uselpg.scst_glm, variables = c("scst", "wave"))
margins_usegrid <- margins(model.dd.vfe.usegrid.scst_glm, variables = c("scst", "wave"))

summary_uselpg <- summary(margins_uselpg) %>% 
  mutate(factor = case_when(factor == "scst" ~ "SC/ST",
                            factor == "wave" ~ "2018"))
summary_usegrid <- summary(margins_usegrid) %>% 
  mutate(factor = case_when(factor == "scst" ~ "SC/ST",
                            factor == "wave" ~ "2018"))

stargazer(summary_uselpg, summary = FALSE, float = FALSE, rownames = F,
          out = "./Manuscript/Tables/uselpg.probit.ame.tex")

stargazer(summary_usegrid, summary = FALSE, float = FALSE, rownames = F,
          out = "./Manuscript/Tables/usegrid.probit.ame.tex")

